Cluster and reentrant anomalies of nearly Gaussian core particles 
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We study through integral equation theory and numerical simulations the structure and dy- 
namics of fluids composed of ultrasoft, nearly Gaussian particles. Namely, we explore the fluid 
phase diagram of a model in which particles interact via the generalized exponential potential 
u(r) = eexp[— (r/a) n ], with a softness exponent n slightly larger than 2. In addition to the well- 
known anomaly associated to reentrant melting, the structure and dynamics of the fluid display 
two additional anomalies, which are visible in the isothermal variation of the structure factor and 
diffusivity. These features are correlated to the appearance of dimers in the fluid phase and to the 
subsequent modification of the cluster structure upon compression. We corroborate these results 
through an analysis of the local minima of the potential energy surface, in which clusters appear 
as much tighter conglomerates of particles. We find that reentrant melting and clustering coexist 
for softness exponents ranging from 2 + up to values relevant for the description of amphiphilic 
dendrimers, i.e., n = 3. . 
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I. INTRODUCTION 

Hard colloidal suspensions have been extensively stud- 
ied as the simplest model of particulate systems. In con- 
trast, several colloidal suspensions are composed of soft, 
macroscopic particles, and display a much more complex 
behavior 0, [2| ■ Star polymers 0-[5( and microgels [Q-Q 
are notable examples of purely repulsive, soft colloids, in 
which softness can be tuned by either altering the macro- 
molecular architecture or by exploiting the sensitivity to 
temperature and solvent properties. As a consequence 
these macromolecular aggregates can be compressed to 
attain very large packin g fr actions — significantly above 
random close packing [i^Tl0| . 

In theory and simulation, soft colloids can be studied 
by modeling them as point particles interacting through 
a soft repulsive potential u(r). Paradigmatic examples 
are Hertzian particles (Tlj . whose repulsion originates 
from the elastic deformation of spheres at contact, and 
the Gaussian core model (GCM) [H, . The thermo- 
dynamic and dynamic behaviors of fluids composed of 
such soft particles differ significantly from the usual hard 
sphere model and display an anomalous, reentrant den- 
sity dependence [llTl25j . Upon increasing density from 
the dilute regime, the particles first pack more densely 
and the structural order increases — as in hard spheres. 
Upon further compression, particles tend to decorrelate 
spatially, since the energetic gain of avoiding overlaps is 
overwhelmed by the increase of available configurations 
(entropic effect) This leads to looser spatial cor- 

relations and enhanced diffusivity From a physical point 
of view, the anomalous behavior can thus be understood 
in terms of the competition between entropic and ener- 
getic effects. The anomalous behavior of the fluid has 
a direct counterpart in the thermodynamic phase dia- 
gram, which displays reentrant melting upon increasing 



density [TJ EH 113 • Very recently, such anomalies have 
indeed been observed in experiments on microgels, which 
can be modeled rather well using Hertzian particles [2?j| . 

On the other hand, a different type of anomaly can 
occur if the pair potential is bounded at the origin, i.e., 
u(r = 0) is finite, but has a steep hill at some length 
scale a. An ideal example is the penetrable sphere model 
(PSM) in which u(r) is a finite constant at r < a and zero 
otherwise [I?} ■ At sufficiently high density, such ultrasoft 
particles display cluster formation. In the high-density, 
fluid phase, the radial distribution function shows a clear 
peak at r = 0, meaning that the particles can interpene- 
trate and form relatively stable clusters p8j . Depending 
on the details of the interaction potential, formation of 
cluster crystals [28U32| and cluster glasses [33| have also 
been observed at high density. In such cluster phases, 
the collective structure is essentially frozen, but indi- 
vidual particles may diffuse by hopping on the cluster 
site s 134, [3 5l| . Through a mean-field analysis, Likos et 
ai. (aalsj have proposed that clustering occurs when- 
ever the Fourier transform of potential, u(k), has some 
negative components (Q^ potentials). 

Both the reentrant anomalies and clustering can be 
studied within a coherent model system, called the gen- 
eralized exponential model (GEM) 



u(r) = e exp[— (r/cr)' 



(1) 
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where n is a softness exponent [29(. Clearly, the GCM 
and the PSM are recovered when n = 2 and n = 00, 
respectively. The GEM has been extensively studied for 
n = 4 and 8, but very little is known about the behav- 
ior for other values of n. According to the mean-field 
argument of Likos et al. [H, H3; clustering occurs for 
n > 2. Thus, the model allows to interpolate continu- 
ously between reentrant anomalies and cluster formation, 
with an interesting crossover expected around n = 2. 
Moreover, values of n around 3 are of interest, as they 
may be used to describe the effective interactions be- 
tween macromolecules with low fractal dimensions, such 
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as amphiphilic dendrimers [38(. Although these effective 
interactions are derived in the infinite dilute regime, a 
recent study indicates that their phenomenology may be 
observed in fully atomistic simulation of dendrimers at 
finite concentrations [39| 

In this paper we present a theoretical and numerical 
investigation of the GEM-n with varying softness n > 2. 
We show that the structural and dynamic properties dis- 
play both reentrant and cluster anomalies when n is 
larger than but close to 2. First, we conduct a prelimi- 
nary investigation of the fluid phase diagram by means 
of integral equation theory for various values of n. We 
focus then on the behavior of the model with n = 2.2. 
We find that the fluid structure and dynamics display 
a sequence of anomalies, which manifests themselves as 
multiple minima and maxima in the isothermal proper- 
ties. We find that the origin of the anomalies observed at 
high density are due to formation of clusters in the fluid. 
To provide a sharp characterization of cluster formation 
we analyze the potential energy surface of the model, fol- 
lowing numerical techniques well-established in the con- 
text of supercooled liquids An analysis of the local 
minima sampled above the melting temperature allows to 
obtain a coherent picture of the structural and dynamic 
anomalies visible at high density. 



II. METHODS 

We consider a model of point particles interacting via 
the generalized exponential potential of order n, Eq. (TTJ . 
As usual, e, <r, and \J ma 2 /e provide the units of length, 
energy, and time, respectively. We use integral equation 
theory to predict the static properties over a broad range 
of state parameters and n values. The Ornstcin-Zernikc 
equation with the hypernettcd chain closure (HNC) re- 
lation [4l| was solved by an iterative algorithm. We 
also used the mean spherical approximation (MSA) clo- 
sure relation [l4|, , which can be solved analytically. 
Molecular dynamics simulations were carried for systems 
composed of N = 4000 particles enclosed in a cubic box 
with periodic boundary conditions. Fluid samples were 
prepared starting from a random configuration using a 
Berendsen thermostat to adjust the temperature during 
the equilibration phase. At low temperature, additional 
runs were carried out using a slow annealing to make sure 
that static and dynamic properties were reproducible and 
were not affected by partial crystallization or phase sep- 
aration. State points for which at least one sample crys- 
tallized were not considered in our analysis, which only 
focused on the fluid portion of the phase diagram. Pro- 
duction runs were performed in the microcanonical en- 
semble using the velocity- Verlet algorithm. The integra- 
tion time step was 0.05 in reduced units for both equili- 
bration and production runs. The properties of the GCM 
are known to be sensitive to the cut off radius r c of the 
potential [111 |42j ■ Therefore, a relatively long range cut 
off was applied (r c = 4.0). In addition, we employed a 



cubic spline interpolation [43| between r c — 4.0 to 4.17 to 
improve energy conservation during the MD simulations 
and to avoid numerical difficulties in the analysis of the 
potential energy surface. The local minima (or inherent 
structures [Hj]) of the potential energy surface were lo- 
cated by minimizing the total potential energy U({ri}) 
starting from independent, instantaneous fluid configu- 
rations. Thanks to the large system size, a relatively 
small number of independent configurations (20-40) was 
sufficient to obtain good quality data. Some larger data 
sets (up to 300 configurations) were considered at spe- 
cific state points. The LBFGS algorithm [45[ was used to 
locate the local minima. We found that convergence be- 
comes particularly slow close to the clustering crossover 
p c (see Sec. IIIip . Depending on density and temperature, 
from 1000 to 50000 iterations were necessary to reduce 
the norm of the total force per particle W ~ -h ff to 
values smaller than 10~ 13 . Such a tolerance criterion was 
sufficient to reduce to zero the number of unstable modes 
of the dynamical matrix. 



III. RESULTS 

The aim of this work is to investigate the interplay be- 
tween reentrant anomalies and clustering in nearly Gaus- 
sian particles, i.e. for values of n slightly larger than 2. 
To understand how the phase diagram is modified as the 
interactions deviate from the pure Gaussian potential, 
we start with a preliminary exploration of fluid structure 
of the model by varying systematically the exponent n. 
For this purpose, we employ the integral equation theory 
using the HNC closure relation, which provides a good 
description of the structure of soft particles at sufficiently 
high density 0,13. 

A simple way to highlight the anomalous features of 
the structure in soft core particles is to identify the locus 
of state points in the T — p plane for which the height 
S* = S(k*) of the first peak of the structure factor is 
constant. We found that HNC suffers from convergence 
problems at intermediate densities, in particular close 
to the clustering crossover (see below), when n < 2.05. 
Therefore, we restrict ourselves to a relatively small value 
of S* . For the models at hand, we use the condition 
S* = 2.3, which is smaller than typical values observed 
at melting in simple li quid s (S* ~ 2.8 according to the 
Hansen- Verlet criterion [43]). We emphasize that the aim 
of this preliminary analysis is to highlight the presence 
of a structural anomaly in the (T, p) diagram and not to 
track closely the fluid-solid phase boundary. Such an op- 
eration would require adjusting the value of S* depend- 
ing on density, since the first peak of S(k) reaches large 
values (i.e., S* > 4.0) along the cluster crystallization 
line |U|37|]- 

The iso-S* curves obtained within the HNC approxi- 
mation are shown in Fig. [T]for several values of n > 2, 
along with the fluid-solid phase boundary of the GCM 
determined in Rcfs. [l|| and [24J. At low density, HNC 
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FIG. 1. Integral equation theory predictions for the iso-S* 
lines in the T-p plane for various values of n. The lines connect 
state points for which the height, S*, of the first peak of 
the structure factor is equal to 2.3. Full and dashed lines 
correspond to HNC and MSA approximations, respectively. 
Open squares indicate the fluid-solid phase boundary of the 
GCM (n = 2) as determined in Refs. Qj] and [13|. For clarity, 
the iso-S 1 * lines predicted by MSA are not shown below the 
fluid-solid phase boundary. 



predicts the iso-S"* line to depend very little on n and to 
follow the low density branch of the GCM phase bound- 
ary. In this portion of the phase diagram, the HNC pre- 
dictions are semi-quantitative at best, since we expect 
the Hansen- Verlet criterion to work well along the low 
density phase boundary. Thus, larger values of S* would 
be expected. For all the values of n shown in the figure, 
the iso-S 1 * lines show a first decrease upon increasing den- 
sity, following the non-monotonic, anomalous behavior of 
the melting line of the GCM. Upon increasing the density 
even further, an additional branch appears: for values of 
n strictly larger than 2, the iso-S* lines follow a common 
envelope until they reach a crossover density and start 
bending upwards. The crossover density diverges when 
n — > 2 + , as it can be inferred by tracking the location 
of the minimum of S* (p) as a function of n (not shown) . 
Such a behavior is a striking confirmation of the mean- 
field argument of Likos et al. (36|. In fact, as it will be 
clear from the following, the minimum of the iso-S* lines 
marks the onset of clustering in the fluid phase. We find 
that the two coexisting anomalies arc visible up to values 
of n around 3, at which the minimum transforms into a 
small inflection. Above n = 3 the iso-S* lines predicted 
by HNC display a purely monotonic increase with den- 
sity. 

Results obtained using the MSA are also included in 
the figure. In this case, the iso-S* lines can be obtained 
analytically and are given by straight lines in the T-p 
plane 



FIG. 2. Fluid phase diagram of the GEM-2.2 in the T-p plane. 
State points at which S* = 2.8 are indicated as full squared 
and are connected by full lines (iso-S 1 * lines). State points at 
which g(0) = 0.05 are drawn as filled or open circles according 
to whether they are characterized by clustering or transient 
overlaps, respectively, as explained in the text. Open squares 
indicate the fluid-solid phase boundary of the GCM (n = 2) 
as determined in Refs. and [24[ . The light grey and dark 
grey shaded areas indicate portions of the phase diagram dis- 
playing finite fractions of dimers and trimers, respectively, in 
the local minima of the potential energy surface. Inset: mag- 
nification of the behavior of the iso-S* line at large densities 
for S* = 2.8 (squares) and S* = 2.3 (triangles). Also included 
are the corresponding results from HNC (solid lines). 



where v* is the value of the minimum of the Fourier trans- 
form of the potential. We see that the MSA predictions 
follow closely those obtained within HNC at high density. 
However, the crossover between the low and intermedi- 
ate density regimes (reentrant anomaly) is not captured 
within the mean-field description of MSA, as the latter 
is only reliable at high temperatures and densities. 

For values of the softness exponents n between 2 + 
and approximately 3, HNC thus predicts the existence 
of three regimes in the fluid phase diagram separated by 
two anomalies: (i) the "reentrant" anomaly well known 
from studies of the GCM and Hertzian spheres, corre- 
sponding to a maximum in the iso-S* lines, and (ii) a 
"cluster" anomaly, observed as the particles start to ex- 
plore the ultrasoft portion of the potential, which appears 
as a minimum in the iso-S* lines. To test this anomalous 
scenario in-depth, we focus in the following on the case 
n = 2.2, for which we collected the most extensive set of 
simulation data. Simulations for different values of n will 
be briefly discussed at the end of this section. 

In Fig. [5] we combine our molecular dynamics results 
and the analysis of the potential energy surface to char- 
acterize the fluid phase diagram the GEM-2.2. We start 
by discussing the shape of the iso-S* line, which we de- 
termined using the value S* = 2.8 (Hansen- Verlet crite- 
rion) . At low density, the line first follows the trend of the 
GCM fluid-solid phase boundary. Upon further comprcs- 
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FIG. 3. Comparison of the structure factors S(k) [panels (a) 
and (b)] and radial distribution functions g(r) [panel (c) and 
(d)] obtained using the HNC closure (solid lines) and MD 
simulations (crosses). Data are shown for p = 0.8, T = 0.01 
in (a) and (c), and for p = 1.4, T = 0.018 in (b) and (d). 



sion, it first passes through a maximum located around 
p r ~ 0.3 and then a minimum around p c m 0.78. Thus 
both the reentrant and cluster anomalies predicted by 
HNC are visible in the simulation data. Upon closer in- 
spection, however, it is clear that the agreement between 
HNC and the simulations in this portion of the phase di- 
agram is only qualitative. The minimum separating the 
soft-core and the cluster regimes is in fact deeper than 
predicted by HNC. These discrepancies are further high- 
lighted in Figs.[J3a) and (c), where we compare the HNC 
predictions for the structure factor S(k) and the radial 
distribution function g(r) around the crossover. HNC 
predicts clustering to occur at smaller densities (or higher 
temperatures) than actually observed in the simulations, 
as seen from the differences in g(r) around r = 0. Only 
at larger densities and temperatures [see Figs. EJb) and 
(d)], the HNC predictions agree quantitatively with the 
simulations, as expected from previous work [3(j. This 
is further illustrated by the 180-5* lines shown in the in- 
set of Fig. [21 where we observe good agreement between 
HNC and simulations for S* = 2.3, but not for S* = 2.8. 
We conclude that better closure relations are required 
for a quantitative description of the fluid phase diagram 
close to the fluid-solid phase boundary, especially at low 
and intermediate densities. 

To illustrate the connection between the minimum of 
the iso-S 1 * line and clustering, we evaluate the radial dis- 
tribution function g(r) and analyze its isothermal varia- 
tion. In the following, we will denote with r* the posi- 
tion of the peak of g(r) corresponding to the first neigh- 
bor (FN) shell. In Fig. HJa) we show the g(r) along 
the isotherm T = 0.012. As the density is increased 
from p = 0.4, the first peak initially shifts markedly to- 
wards smaller distances, as expected for soft-core parti- 
cles. Upon crossing p c , a broad peak develops at r = 



and grows by further increasing the density. The pres- 
ence of a well-defined minimum between r = and r* , 
indicates the formation relatively stable clusters. The 
FN peak position r* also shows an interesting behavior: 
it first decreases with increasing density up to p c , then it 
starts to increase until the density reaches approximately 
p = 1.1. This point will be further discussed below. 

In contrast to cluster formation, transient overlaps can 
always occur at any density in fluids of ultrasoft particles 
as a consequence of thermal fluctuations. Indeed, at high 
temperatures, coreless particles often superpose during 
collisions, which implies a finite value of g(r = 0). How 
to distinguish between clustering and such transient over- 
laps? We propose to classify clustering states as those for 
which the second derivative g"(r = 0) < 0. In practice, 
an equivalent and numerically more stable criterion is to 
check whether the g(r) displays a minimum between 
and r*. States for which g(0) is finite but <?"(0) > arc 
characterized instead by transient overlaps. 

With all this information at hand we return to the 
fluid phase diagram Fig. [5] and mark the boundary be- 
tween normal and cluster states. To this end, we identify 
the state points in the T-p plane at which 17(0) reaches 
constant, "small" value. The reference value g(0) = 0.05 
is small enough to detect the early onset of clusters and 
overlaps, still sufficiently large to minimize statistical un- 
certainties. Numerically, we evaluate g(0) from the prob- 
ability of finding particles at distances less than 0.03. Fol- 
lowing the classification introduced above, we use open 
or filled circles in Fig. [2] for states characterized by clus- 
ters or transient overlaps, respectively. In the low and 
intermediate density regime, the onset of transient over- 
laps is located at temperatures well above the melting 
line. Upon increasing the density, the iso-g(O) line bends 
downwards and starts to run vertically in the T — p plane. 
In this portion of the phase diagram, the sign of g"(0) 
along iso-g(O) line becomes negative, signaling the tran- 
sition to cluster states. By visual inspection, it is clear 
that the boundary between normal and cluster states cor- 
responds well to the minimum of the iso-5 1 * line, which 
motivates the name "cluster anomaly" . 

To understand the physical origin of the cluster 
anomaly, we analyze in more detail the changes in the 
microscopic structure of the clusters. To identify clus- 
ters, we use the same geometrical definition as in pre- 
vious work [33|, [3f| : a particle is considered to be part 
of a cluster if its distance to at least another particle of 
that cluster is smaller than a fixed value d. Typically, 
d is chosen around the first minimum of g(r), which is 
located in the range 0.6-0.7 over the pertinent range of 
density and temperatures. For simplicity, in our analysis 
we used a state independent cut off d = 0.6. 

The cluster population = (Nk)/N, where Nf- is the 
number of fc-mers in a configuration, is shown in Fig.^b) 
as a function of density along the isotherm T = 0.012. 
We find that starting from p w 0.78 the fluid is no longer 
composed only of isolated particles (monomers) but pos- 
sesses a growing fraction of dimers. In the range of den- 
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FIG. 4. Density variation of the radial distribution function and cluster properties of the GEM-2.2 evaluated for instantaneous 
fluid configurations: (a) g(r) along the isotherm T — 0.012. For clarity, each line has been vertically shifted from the previous 
one by a constant offset, 1.25. For p > 0.8, the horizontal dotted lines indicate the ordinates at which g(r) would equal zero. 
Filled circles indicate the first neighbor shell peaks, (b) Fraction of fc-mers for T = 0.012 as a function of p for k — 1 (monomers, 
triangles), 2 (dimers, circles), and 3 (trimers, squares), (c) and (d) Monomer- monomer (dotted line), monomer-dimer (dashed 
line) and dimer-dimer (solid line) radial distribution functions evaluated at (c) p = 0.8, T = 0.012 and (d) p = 0.9, T — 0.012. 
Note that the distribution functions are zero by construction for r < d = 0.6. 



sity between 0.78 and 1.1, we observe coexistence be- 
tween monomers and dimers, the fraction of higher order 
fc-mers being zero within our numerical accuracy. We re- 
mark that this range varies slightly as a function of the 
chosen temperature. To illustrate the real space structure 
of coexisting monomers and dimers, we show in Fig|4jc) 
and (d) the monomer- monomer (MM), monomer-dimer 
(MD) and dimer-dimer (DD) radial distribution func- 
tions evaluated at p = 0.8 and p = 0.9, respectively. 
These correlation functions are calculated from 

y IN a N P V 

where a and /3 are either M or D and the term i = j 
in the double sum is ignored when a = [3. For dimers, 
the distance fij is evaluated from their center of mass. 
We find that monomers are more loosely correlated than 
dimers, as it can be evinced from the nearly flat shape 
of (?mm (r). This probably reflects the fact the effective 
interactions between dimers are harsher than the bare 
potential between individual particles. Dimers, on the 
other hand, tend to occupy a larger volume and are thus 
separated from their neighbors by larger distances. This 
suggests that the slight increase of r* above p c is due to 
the larger FN shell of the growing dimers population. 

Upon further compression above p = 1.1, a small frac- 
tion of trimers P3 appears and grows by increasing den- 
sity at the expenses of P\ and P2 [see Fig. Efb)]. Keep- 



ing T — 0.012 fixed, we cannot extend the analysis of 
the fluid regime beyond p ~ 1.3, due to incipient, par- 
tial cluster crystallization. We thus conclude that, above 
p c , clusters are characterized by low population numbers, 
i.e. 2 or 3, and coexist in a nearly structureless contin- 
uum formed by isolated particles. This scenario is most 
likely connected to the underlying crystalline phase di- 
agram, which may possess a cascade of transitions in- 
volving normal fluids (i.e., isolated particles) and cluster 
crystals with integer population numbers. Such a cas- 
cade of complex crystal structures is known to occur in 
the PSM Hi, [32| and in the GEM-4 HHH]. 

Due to the thermal motion of particles, however, the 
definition of a cluster in the fluid states is subject to some 
ambiguity. For instance, the drop of <?mm(?") at r = d for 
p = 0.9 (and to a less extent for p = 0.8) indicates that 
there is no clear separation between the monomers and 
the FN shells of the clusters. To obtain a sharper picture, 
we apply our cluster analysis to the local minima of the 
underlying potential energy surface (PES). We perform 
a statistical sampling of the PES properties by quench- 
ing several independent samples (instantaneous configu- 
rations) at a given p and T, as described in the Sec.|Tl] In 
the following, we will focus on the isotherm T = 0.02 and 
study how the underlying landscape changes upon vary- 
ing the density. Although the PES is fully determined 
once the number of particles and the volume are specified, 
the properties of the basins visited by the system may 
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FIG. 5. Density variation of the radial distribution function and cluster properties of the GEM-2.2 evaluated for local minima 
of the potential energy surface. Lines and symbols as in Fig. [4] Note that in this case both the radial distribution function 
j mm (r) and the cluster populations are independent of T over the studied range of state parameters. The actual data shown 
in the figure were obtained from minimizations at T — 0.02. In (a), solid lines are colored according to the different regimes 
of cluster populations, also indicated by the shaded areas in (b) . Vertical arrows highlight the splitting of the FN maxima at 
p — 0.75 and p — 1.3 and indicate the location of the secondary peaks. In (c) and (d) the partial radial distribution functions 
are shown for p = 0.8, T = 0.02 and p = 0.9, T = 0.02, respectively. 



vary below some onset temperature [40j . However, this 
is not observed for the GEM-2.2, at least down to tem- 
peratures close to the fluid-solid phase boundary. The 
situation resembles the scenario observed in simple liq- 
uids, where the properties of local minima visited above 
melting are indeed temperature independent [40l |. 

The main advantage over the analysis based on instan- 
taneous configurations is that, in the local minima, clus- 
ters appear as much tighter conglomerates of particles, 
sitting practically on top of each other. This is demon- 
strated by the radial distribution function g min (r) eval- 
uated using the minimized configurations, see Fig. [5ja). 
The peak at zero distance, which emerges at the cluster- 
ing crossover, is now extremely narrow and is separated 
from the FN shell by a clear gap where g(r) is zero within 
numerical precision. This, in turn, makes the definition 
of clusters straightforward and unambiguous. In prac- 
tice, we used d = 0.2 as a cut off distance for the cluster 
analysis on local minima. 

Compared to instantaneous configurations, local min- 
ima present a more subtle density dependence of the ra- 
dial distribution function. The latter is best understood 
when discussed along with that of the cluster population 
P™ m , which we show in Fig. [SJb) . An analysis of the clus- 
ter populations reveals a steep rise of dimers at p c = 0.78, 
making the definition of the clustering crossover very 
sharp. The fraction of trimers becomes non zero only 
above p = 1.1, again displaying a sharper growth than in 



instantaneous configurations. Upon increasing the den- 
sity even further (p > 1.3), monomers eventually dis- 
appear and we observe coexistence between dimers and 
trimers. These various regimes are indicated in Fig. [2] 
by shaded areas, to highlight their connection to the ob- 
served anomalies. 

We can now more easily understand the evolution of 
g min (r) with density. The FN peak of g min (r) shows a 
broadening around p c and then splits upon further in- 
creasing the density. The two peaks arise from the coex- 
istence of monomers and dimers, which are characterized 
by two distinct typical inter-cluster separations. This is 
confirmed by an analysis of the monomer-monomer and 
dimer-dimer contributions to g min (r), which are calcu- 
lated as for instantaneous configurations and are shown 
in Fig. [SJc) and (d) at p = 0.8 and p = 0.9, respec- 
tively. In particular, the splitting of the FN shell for 
p = 0.75 [see panel (a)] is due to growing dimer-dimer 
correlations. Indeed, the position of the secondary max- 
imum at r rs 1.45 matches the location of the first peak 
in <7dd (r). We also find that the MM peak at small dis- 
tances soon looses its structure above p c , whereas the 
DD peak shifts progressively to smaller distances. This, 
together with the density evolution of g min (r) above p c , 
suggests that dimers behave as larger, effective soft-core 
particles. Above p = 1.1, the FN peak of g nun (r) shows 
again a broadening and a subsequent splitting, which 
arises now from the coexistence of subpopulations of 
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dimcrs and trimers in the local minima. 

The partial radial distribution functions g™™(V) for 
p = 0.9 [see Fig. Efd)] qualitatively confirm the re- 
sults obtained from instantaneous configurations in the 
monomcr-dimcr coexistence region. In particular, we see 
that isolated particles, which amount to ~ 10% of the 
fc-mers at this density, are spatially uncorrelated, as in- 
dicated by the unstructured profile of gwwi{r). Very close 
to the cluster crossover [see Fig.[5jc)] , we observe a strong 
suppression of g^§{r) and and a concomitant enhance- 
ment of monomer- monomer correlations. This indicates 
that monomers and dimers have a tendency to microseg- 
regate and suggests an underlying instability with re- 
spect to phase separation into coexisting fluid and clus- 
ter phases. Visual inspection of the real space structure 
of the local minima at p = 0.8 confirms these observa- 
tions. We found no such tendency to microsegregation in 
the instantaneous configurations. These results resemble 
the scenario observed in certain phase separating binary 
Lennard- Jones mixtures [49| , for which the local minima 
display clear signs of demixing — even for states at which 
the equilibrium fluid is essentially homogeneous. 

We now analyze the connection between the structural 
anomalies and the dynamics. It is well established [l9[ 
that the reentrant anomaly in soft-core models is accom- 
panied by a dynamic anomaly in the single particle dif- 
fusivity: beyond the reentrant crossover, the dynamics 
accelerates upon increasing density due to the progres- 
sive the loss of structural correlations. We now show that 
also the cluster anomalies also have their dynamic coun- 
terparts. In the three panels of Fig. [5] we compare the 
isothermal variation of the peak height S* of the struc- 
ture factor, the inverse of the two-body contribution to 
the excess entropy 

/>OG 

= -2tt P / r 2 {g(r) \ng(r) - [g(r) - 1]} dr 
Jo 

and the diffusion constant D, which has been obtained 
from the usual Einstein relation 
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FIG. 6. Structural and dynamic anomalies of the GEM-2.2 
as a function of density along selected isotherms: T = 0.02 
(squares), 0.014 (circles), 0.012 (triangles), 0.011 (reversed 
triangles), and 0.01 (diamonds), (a) Height, S* , of the first 
peak of the structure factor, (b) Inverse of two-body excess 
entropy s™. (c) Diffusion constant D. Arrows at p — 0.78 
and p = 1.1 mark the appearance of dimers and trimers in 
the PES, respectively. 



(R 2 (t)) -> 6Dt 

where R 2 (t) = Y^i=\ 1^(0 — ^(0)| 2 ' 1S the mean square 
displacement. The two-body excess entropy s| x is a sim- 
ple measure of pair translation order in the system. This 
quantity has been used to characterize structural anoma- 
lies in fluids of ultrasoft particles 0, H(| ■ 

The two isotherms at high temperature reproduce the 
reentrant anomalies in both structure and dynamics, as 
expected. On the contrary the two low temperature 
isotherms necessarily hit, at small p, on the fluid-bec 
phase boundary, and therefore only the cluster anomalies 
are visible. We find that s| x follows closely the density 
variation of the S* . A close inspection of panels (a) and 
(b) of Fig. [6] reveals only a minor discrepancy between 
the location of the minima in S* and s| x , with the latter 
providing a slightly larger 5%) value of the density 
at the minimum. A comparison of the structural metrics 



and D shows that the clustering crossover at p c is clearly 
accompanied by a decrease in diffusivity. Physically, this 
anomaly can be easily understood, since transient clus- 
ters temporarily behave as effective, heavier particles, 
thus slowing down the dynamics of the fluid. Despite 
the clear connection between structure and dynamics il- 
lustrated by Fig. |H1 an attempt to collapse the dynamic 
data as a function of the two-body excess entropy sf 1 , as 
in previous works [III [2(| , did not yield satisfactory re- 
sults on a quantitative level (see also [50j | for a discussion 
on the breakdown of such a "Rosenfeld scaling" j5lj). 

Around p w 1.1 and at sufficiently low temperatures, 
an additional anomaly is visible as a maximum in S* 
(and s| x ) and a minimum in D. The relatively large val- 
ues attained by S* along the lowest isotherm, T = 0.01, 
suggest that the system might be slightly supercooled, 
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FIG. 7. Structural and dynamic anomalies of the GEM-3 
as a function of density along selected isotherms: T = 0.07 
(squares), 0.05 (circles), 0.04 (triangles), 0.03 (reversed trian- 
gles), and 0.02 (diamonds), (a) S*, of the first peak of the 
structure factor, (b) Inverse of the two-body excess entropy 
s™. (c) Diffusion constant D. 



despite being (meta) stable during the whole length of 
our simulations. Figure Uta) and (b) show however that 
this maximum develops smoothly by decreasing temper- 
ature. The inflection of the iso-5 1 * lines in the (T, p) 
diagram [see inset of Fig. [5] is most likely a high temper- 
ature manifestation of this additional anomaly. Thus, 
we are confident that this represents an actual feature 
of the equilibrium fluid. This anomalous scenario is also 
observed at smaller n values, down to at least 2.01 (not 
shown). In this case, the two cluster anomalies are visi- 
ble at much lower temperatures than in the GEM- 2. 2, as 
expected from the HNC results. 

Qualitatively, the second cluster anomaly in the GEM- 
2.2 correlates with the appearance of trimers and the 
disappearance of isolated particles. However, its under- 
lying physical origin remains unclear: for instance, why 
should the formation of trimers accelerate the dynam- 
ics? A possible explanation is that, at these densities, 
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FIG. 8. Fraction of fc-mers in local minima of the potential 
energy surface as a function of n — 2 for k — 1 (triangles), 2 
(circles), 3 (squares), and 4 (reversed triangles). Local min- 
ima were sampled at p = 2.0, T = 0.1 . 



different subpopulations of fc-mers coexist (fc = 1,2,3), 
which enhances the structural disorder and thus acceler- 
ates the dynamics. It could also be argued that dimcrs 
behave as effective, ultrasoft particles and display an ad- 
ditional reentrant anomaly in their own right. We note, 
in passing, that neither HNC nor MSA can reproduce this 
anomalous feature. Further work is definitely needed to 
clarify this issue. 

By varying the softness exponent, we found that co- 
existence of reentrant and cluster anomalies is well vis- 
ible up to at least n = 3, a value that can be used 
to model the effective interactions between amphiphilic 
dendrimers [38[ . This is demonstrated in Fig. [7J where 
we show the isothermal variation of S* , s| x and D in 
the GEM-3. In this case, however, only the first cluster 
anomaly is visible in the range of densities and temper- 
atures covered by our simulations. The absence of the 
second anomaly is consistent with the fact the fraction 
of trimers in the fluid is negligible (P3 < 10 4 ) at the state 
points shown in Fig. [7J 

Finally, we study the evolution of the cluster popula- 
tion in the local minima as n — > 2 + . Figure [8] shows the 
fraction of fc-mers Pk with k = 1, 2, 3, 4 at fixed p = 2.0 
and T = 0.1. As the interaction smoothly transforms into 
a Gaussian, clusters of high order disappear progressively. 
At the studied density, only models with n — 2 < 0.003 
are fully dissociated within our statistics. The variation 
of Pk with n thus provides a striking confirmation of the 
Likos criterion for clustering in fluid states. 



IV. CONCLUSIONS 

In this work we demonstrated the existence of mul- 
tiple anomalies in the structure and dynamics of fluids 
composed of nearly Gaussian particles. We focused in 
particular on particles interacting via the generalized ex- 
ponential potential, Eq. ([T]) , with an exponent n slightly 
larger than 2. Our most extensive set of data concerns 
the GEM-2.2, for which we gathered molecular dynam- 
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ics results and for which we analyzed systematically the 
potential energy surface. 

At low densities, the structure and dynamics of the 
fluid display a non-monotonic density variation. Such 
a reentrant anomaly is well-known in soft-core parti- 
cles [HI EH and manifests itself as a loss of translation 
order upon compression and a concomitant acceleration 
of the single-particle dynamics. Nearly Gaussian parti- 
cles show two additional anomalies by compressing the 
fluid at constant temperature: the height S* of the first 
peak of the structure factor displays a minimum and an 
additional inflection at higher densities. At sufficiently 
low temperatures, such an inflection transforms into a 
small but well-defined minimum. All these anomalies 
have counterparts in the single particle dynamics, with 
minima of S* corresponding to maxima of the diffusion 
constant D, and viceversa. We note that a simple scaling 
between Newtonian and overdamped Langevin dynamics 
has been recently reported for ultrasoft fluids [521 ]. indi- 
cating that the connection between structural and dy- 
namical anomalies observed here is likely quite general. 

At high densities, particles explore the bounded part 
of the potential. The two additional anomalies occur- 
ring in this portion of the fluid phase diagram, can thus 
be attributed to the appearance of clusters in the fluid. 
Namely, we found that the minimum in S* (p) correlates 
with the sharp growth of dimers in both instantaneous 
configuration and in the local minima of the PES. The 
additional maximum ostensibly correlates to the appear- 
ance of trimers upon further compression, although its 
actual physical origin needs to be clarified. Interestingly, 
the spatial structure of the local minima reveals that co- 
existing subpopulations of isolated particles and clusters 
can show a slight segregation. This suggests an underly- 
ing instability of fluid with respect to phase separation 
into coexisting fluid and cluster phases. It would thus be 
interesting to characterize in more detail the low temper- 
ature phase diagram of nearly Gaussian particles, along 
the lines of recent numerical work on the GEM-4 [3l[ . 

We demonstrated that clustering in the fluid phase dis- 
appears in a continuous way as the interaction potential 
tends to a purely Gaussian shape, i.e., as the softness 
exponent n approaches 2 from above. This is an ex- 
plicit confirmation of the mean-field argument by Likos 
ct al. (36[ , which has never been so far tested explicitly 



on GEMs with n close to 2. What our results clarify, 
however, is that reentrant and clustering behavior are 
not mutually exclusive, as is sometimes implicit in the 
presentation of Q + and models (36|. We found that 
coexistence of reentrant and cluster anomalies is well vis- 
ible up to at least n = 3. This behavior was anticipated 
by Zhang et al [31| in their study of the low tempera- 
ture phase diagram of the GEM-4 model and has been 
confirmed here more systematically. Whether these coex- 
isting anomalies may be observed in colloidal suspensions 
of highly branched macromolecules, such as dendrimers, 
is an open question that awaits an experimental verifica- 
tion. 

Another interesting implication of our work concerns 
the nature of the possible glassy states of ultrasoft par- 
ticles. It has been shown that the GCM becomes a sur- 
prisingly goo d glass-former at sufficiently high densities 
(p ~ 2.0) [23|]. This has been attributed to a frustration 
effect induced by the long range of the Gaussian poten- 
tial. We found that for n = 2.01 the cluster anomaly is 
shifted up to approximately p ~ 2.0. This, in turn, opens 
up a broad region in the fluid phase diagram in which the 
model behaves essentially as the GCM. In such a region, 
models with n values marginally larger than 2 might dis- 
play peculiar glass-forming properties. In particular, it 
would be interesting to see whether the proximity to clus- 
ter formation destabilizes the glassy behavior of the GCM 
and to what extent such a behavior will be affected by 
clustering at even higher densities. Finally, we believe 
that the analysis of the PES may prove particularly use- 
ful in the study of amorphous cluster phases of size dis- 
persed ultrasoft particles [33|. Work to investigate the 
properties of the PES explored by cluster glass-forming 
fluids is currently underway. 
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